其他
临床预测模型之二分类资料ROC曲线绘制
ROC曲线是评价模型的重要工具，曲线下面积AUC可能是大家最常见的模型评价指标之一。
如果你还不太了解关于ROC曲线中的各种指标，请看下面这张图，有你需要的一切（建议保存）：
R语言中有非常多的方法可以实现ROC曲线，但是基本上都是至少需要2列数据，一列是真实结果，另一列是预测值，有了这两列数据，就可以轻松使用各种方法画出ROC曲线并计算AUC。
这篇文章带大家介绍最常见的并且好用的二分类变量的ROC曲线画法。
方法1
方法2
方法3
方法1
使用pROC包，不过使用这个包需要注意，一定要指定direction，否则可能会得出错误的结果。
这个R包计算AUC是基于中位数的，哪一组的中位数大就计算哪一组的AUC，在计算时千万要注意！
使用pROC包的aSAH数据，其中outcome列是结果变量，1代表Good，2代表Poor。
library(pROC)
## Type 'citation("pROC")' for a citation.
##
## 载入程辑包：'pROC'
## The following objects are masked from 'package:stats':
##
## cov, smooth, var
data(aSAH)
dim(aSAH)
## [1] 113 7
str(aSAH)
## 'data.frame': 113 obs. of 7 variables:
## $ gos6 : Ord.factor w/ 5 levels "1"<"2"<"3"<"4"<..: 5 5 5 5 1 1 4 1 5 4 ...
## $ outcome: Factor w/ 2 levels "Good","Poor": 1 1 1 1 2 2 1 2 1 1 ...
## $ gender : Factor w/ 2 levels "Male","Female": 2 2 2 2 2 1 1 1 2 2 ...
## $ age : int 42 37 42 27 42 48 57 41 49 75 ...
## $ wfns : Ord.factor w/ 5 levels "1"<"2"<"3"<"4"<..: 1 1 1 1 3 2 5 4 1 2 ...
## $ s100b : num 0.13 0.14 0.1 0.04 0.13 0.1 0.47 0.16 0.18 0.1 ...
## $ ndka : num 3.01 8.54 8.09 10.42 17.4 ...
计算AUC及可信区间：
res <- roc(aSAH$outcome,aSAH$s100b,ci=T,auc=T)
## Setting levels: control = Good, case = Poor
## Setting direction: controls < cases
res
##
## Call:
## roc.default(response = aSAH$outcome, predictor = aSAH$s100b, auc = T, ci = T)
##
## Data: aSAH$s100b in 72 controls (aSAH$outcome Good) < 41 cases (aSAH$outcome Poor).
## Area under the curve: 0.7314
## 95% CI: 0.6301-0.8326 (DeLong)
plot(res,legacy.axes = TRUE)
可以显示最佳截点，比如AUC最大的点：
plot(res,
legacy.axes = TRUE,
thresholds="best", # AUC最大的点
print.thres="best")
可以显示AUC的可信区间：
rocobj <- plot.roc(aSAH$outcome, aSAH$s100b,
main="Confidence intervals",
percent=TRUE,ci=TRUE,
print.auc=TRUE
)
## Setting levels: control = Good, case = Poor
## Setting direction: controls < cases
ciobj <- ci.se(rocobj,
specificities=seq(0, 100, 5)
)
plot(ciobj, type="shape", col="#1c61b6AA")
plot(ci(rocobj, of="thresholds", thresholds="best"))
多条ROC曲线画在一起：
rocobj1 <- plot.roc(aSAH$outcome, aSAH$s100,percent=TRUE, col="#1c61b6")
## Setting levels: control = Good, case = Poor
## Setting direction: controls < cases
rocobj2 <- lines.roc(aSAH$outcome, aSAH$ndka, percent=TRUE, col="#008600")
## Setting levels: control = Good, case = Poor
## Setting direction: controls < cases
legend("bottomright", legend=c("S100B", "NDKA"), col=c("#1c61b6", "#008600"), lwd=2)
两条ROC曲线的比较，可以添加P值：
rocobj1 <- plot.roc(aSAH$outcome, aSAH$s100,percent=TRUE, col="#1c61b6")
## Setting levels: control = Good, case = Poor
## Setting direction: controls < cases
rocobj2 <- lines.roc(aSAH$outcome, aSAH$ndka, percent=TRUE, col="#008600")
## Setting levels: control = Good, case = Poor
## Setting direction: controls < cases
legend("bottomright", legend=c("S100B", "NDKA"), col=c("#1c61b6", "#008600"), lwd=2)
testobj <- roc.test(rocobj1, rocobj2)
text(50, 50, labels=paste("p-value =", format.pval(testobj$p.value)), adj=c(0, .5))
方法2
使用ROCR，如果你只是为了画一条ROC曲线，这是我最推荐的方法了，美观又简单！
library(ROCR)
使用非常简单，3句代码，其中第2句是关键，可以更改各种参数，然后就可以画出各种不同的图形：
pred <- prediction(aSAH$s100b,aSAH$outcome)
perf <- performance(pred, "tpr","fpr")
auc <- round(performance(pred, "auc")@y.values[[1]],digits = 4)
plot(perf,lwd=2,col="tomato")
abline(0,1,lty=2)
legend("bottomright", legend="AUC of s100b: 0.7314", col="tomato", lwd=2,bty = "n")
添加箱线图：
perf <- performance(pred, "tpr", "fpr")
perf
## A performance instance
## 'False positive rate' vs. 'True positive rate' (alpha: 'Cutoff')
## with 51 data points
plot(perf,
avg="threshold",
spread.estimate="boxplot")
还可以绘制PR曲线，召回率recall为横坐标，精确率precision 为纵坐标：
perf <- performance(pred, "prec", "rec")
plot(perf,
avg= "threshold",
colorize=TRUE,
lwd= 3,
main= "Precision-Recall plot")
plot(perf,
lty=3,
col="grey78",
add=TRUE)
还可以把特异度为横坐标，灵敏度为纵坐标：
perf <- performance(pred, "sens", "spec")
plot(perf,
avg= "threshold",
colorize=TRUE,
lwd= 3,
main="Sensitivity/Specificity plots")
plot(perf,
lty=3,
col="grey78",
add=TRUE)
这个包还可以计算非常多其他的指标，各种图都能画，大家可以自己探索。
方法3
使用tidymodels。这个包很有来头，它是R中专门做机器学习的，我很快就会详细介绍它，它也是目前R语言机器学习领域两大当红辣子鸡之一！另一个是mlr3。
suppressPackageStartupMessages(library(tidymodels))
它很优雅，如果你要计算AUC，那么就是roc_auc()函数：
aSAH %>% roc_auc(outcome, s100b,event_level="second")
## # A tibble: 1 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 roc_auc binary 0.731
如果你是要画ROC曲线，那么就是roc_curve()函数：
aSAH %>% roc_curve(outcome, s100b,event_level="second") %>%
ggplot(aes(x = 1 - specificity, y = sensitivity)) +
geom_path(size=1.2,color="firebrick") +
geom_abline(lty = 3) +
coord_equal() +
theme_bw()
还有太多方法可以画ROC了，不过pROC和ROCR基本上技能解决99%的问题了。
最后，给大家看看cran中比较常见的画ROC曲线的包，大家有兴趣可以自己探索：
library(pkgsearch)
rocPkg <- pkg_search(query="ROC",size=200)
rocPkgShort <- rocPkg %>%
filter(maintainer_name != "ORPHANED") %>%
select(score, package, downloads_last_month) %>%
arrange(desc(downloads_last_month))
head(rocPkgShort,20)
## # A data frame: 20 x 3
## score package downloads_last_month
## * <dbl> <chr> <int>
## 1 4711. caTools 172954
## 2 10868. pROC 123923
## 3 932. ROCR 82311
## 4 1551. plotROC 9719
## 5 2471. PRROC 7209
## 6 282. riskRegression 6622
## 7 173. InformationValue 4710
## 8 1950. survivalROC 4096
## 9 421. PresenceAbsence 3534
## 10 2068. cvAUC 3144
## 11 1594. timeROC 2570
## 12 779. risksetROC 2253
## 13 316. mlr3viz 2243
## 14 109. RcmdrPlugin.EZR 1953
## 15 198. logcondens 1947
## 16 427. ROCit 1744
## 17 1469. precrec 1497
## 18 183. WVPlots 1491
## 19 205. hmeasure 1013
## 20 67.4 wrProteo 965
方法太多了，千万不要迷失哦~
以上就是今天的内容，希望对你有帮助哦！欢迎点赞、在看、关注、转发！
欢迎扫描二维码加 QQ群 ：613637742
欢迎关注公众号：医学和生信笔记
往期回顾
一文搞懂临床预测模型的评价！
简单易懂：什么是临床预测模型？
应用预测建模：过度拟合和模型调优
ggplot2版本的热图-方便拼图！
使用 aplot “拼”一个热图
超详细的R语言热图之complexheatmap系列07